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Abstract 

Recent developments in fracture mechanics analyses of the interfacial crack problem are 
reviewed. The intent of the review is to renew the awareness of the oscillatory singularity 
at the crack tip of a bimaterial interface and the problems that occur when calculating 
mode mixity using numerical methods such as the finite element method in conjunction 
with the virtual crack closure technique. Established approaches to overcome the non- 
convergence issue of the individual mode strain energy release rates are reviewed. In the 
recent literature many attempts to overcome the nonconvergence issue have been 
developed. Among the many approaches found only a few methods hold the promise of 
providing practical solutions. These are the resin interlayer method, the method that 
chooses the crack tip element size greater than the oscillation zone, the crack tip element 
method that is based on plate theory and the crack surface displacement extrapolation 
method. Each of the methods is validated on a very limited set of simple interface crack 
problems. However, their utility for a wide range of interfacial crack problems is yet to 
be established. 


Introduction 

Over the past two decades, the use of fracture mechanics has become common practice 
for characterizing the onset and growth of delaminations [1-3]. In order to predict 
delamination onset or growth, the calculated strain energy release rate components are 
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compared to interlaminar fracture toughness properties which are measured over a range 
from pure mode I loading to pure mode II loading. Figure 1 shows such data (blue circles) 
for a typical carbon/epoxy material. Failure is expected when, for any given mixed mode 
ratio GiiIGt, the calculated total energy release rate, GV, exceeds the interlaminar fracture 
toughness, G c , shown as a curve fit through the experimental data [4], 

Crack onset, propagation, and growth at bimaterial interfaces are of interest to a large 
community. In addition to increasing applications for the analysis composites [5-7], 
cracking has recently been investigated in conjunction with microelectronics [8-11], and 
piezoelectric bimaterials [12] as well as thin film interfaces [13]. Also, cracks that 
develop at the fiber-matrix interface in composites [14], Cracks that propagate 
perpendicular to a bimaterial interface have also been of interest [15] [16]. In the current 
paper the attention, however, is limited to methods that are suitable to predict 
delamination onset, propagation and growth in composites, as well as face sheet/core 
disbonding in sandwich structures. 

The objective of this paper is to renew the awareness of the oscillatory singularity at 
the crack tip at a bimaterial interface and the problems that occur when calculating mode 
mixity using numerical methods such as the finite element method in conjunction with 
the virtual crack closure technique (VCCT) [17-19]. The mathematical background of the 
oscillatory singularity at the crack tip at a bimaterial interface is discussed briefly. 
Established approaches to overcome the non-convergence issue of obtaining the 
individual mode strain energy release rates are reviewed and example cases are 
presented. New methods to provide straightforward practical solutions to easily obtain 
converged solutions for engineering applications are discussed and recommendations for 
future research are presented. 

The paper is organized as follows: First, an overview of the bimaterial interface crack 
problem is given. Second, a brief overview of the virtual crack closure technique together 
with new developments in the methodology is presented. Third, recent developments in 
the analyses of interface crack problems are briefly discussed. Fourth, a more detailed 
description of several methods to obtain satisfactory mode separation is presented. A 
discussion follows and the paper ends with some concluding remarks. 


Bimaterial Interface Crack 

An interface crack of length 2a between two materials A and B is shown in Figure 2, 
subjected to remote normal Sn and tangential tractions St. For example, the two materials 
A and B in Figure 2 may be the neighboring plies in a composite laminate. 


Williams [20] observed that the singularity at the tips of the crack along the interface 
between two dissimilar isotropic materials is of the form Q + iyj where 
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with Kj=3-4Vj is for plane strain conditions, u, is the shear modulus, V/ is the Poisson's 
ratio for j th material and j= A, B. In equation (1), /i is known as the Dundurs parameter 
and is defined in terms of Uj and k-, as 


B = BaO<b ~ 1) -li B (K A - 1) 

Ha(k b + 1) + h b (k a + 1) 

For orthotropic and anisotropic materials the singularity is of similar form and the y 
expression is much more complex than the isotropic case. 


The stresses, t , along the bondline (on the 6 = 0 line) can be written in a form [21-23] 
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and the relative displacements, ur behind the crack tip (at x = a) as 
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where / = . In equations (3) and (4), K is the complex stress-intensity factor defined 

as [21] 


K = (fc x + i fc 2 )V7rcosh(7ry) (5) 

The near-field solution in equations (3) and (4) suggests that the stresses very near the 
crack tip, (r-»0), oscillate (because of the r lY term) and the crack faces interpenetrate 
each other in this region. 


The mode-mixity phase angle, 1 p K , can be expressed as 


1 p K = arctan 


Im (Kr lY ) 
Re(Kr iY ) 


( 6 ) 


where Im and Re are the imaginary and real parts, respectively, of the complex function 
in the brackets. 


Strain Energy Release Rates 

The strain energy release rate, Gt, can be calculated using Irwin’s crack closure integral 
over a virtual crack closure length Aa, from equations (3) and (4), as 
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where K is the complex conjugate of K. The strain energy release rate in equation (7) 
will be denoted as the total strain energy release rate, Gt. While the total strain energy 
release rate appears to control fracture of bimaterial plates, the mode I component of the 
strain energy release rate appears to correlate the experimental data better in fatigue type 
loading [24], Also in composite materials a strong dependence of the experimentally 
determined fracture toughness G c on the individual mode strain energy release rates was 
observed as already shown in Figure 1. Therefore, in the literature considerable attention 
was devoted to the calculation of the individual mode strain energy release rates [18, 19, 
24-27]. Raju et al. [26] showed that for an interface crack between two dissimilar 
isotropic materials, the individual strain energy release rates take the form 

G,= lim Re fC + D • (Aa) ty l 

Aa->0 L J 

G u = lim Re \C - D • (Aa)^l 

Aa->0 L J 

G t = G t + G u = lim Re [2 C] (8) 

Aa->0 

and 

tan 2 x/jq = jf- (9) 

b i 

where Aa is the virtual crack closure length as shown in Figure 3. The parameters C and 
D in equation (8) are complex constants. Equation (8) suggests that the individual strain 
energy release rates, Gy and Gu, depend on Aa and have no well defined limits, while the 
total energy release rate, Gt, is independent of Aa and has a very well defined limit. For 
homogeneous, isotropic materials, however, the results for mode mixity based on 
equations (6) and (9) are identical. 


The Virtual Crack Closure Technique : Overview and Recent Developments 

The virtual crack closure technique [17-19] is widely used for computing energy 
release rates based on results from continuum two-dimensional (2D) and solid three- 
dimensional (3D) finite element (FE) analyses. The mode I and mode II components of 
the strain energy release rate, G/ and Gu respectively are computed using VCCT as 
shown in Figure 4a for a 2D four-node element. For geometrically nonlinear analysis 
where large deformations may occur, both forces and displacements obtained in the 
global coordinate system need to be transformed into a local coordinate system (x', y) 
which originates at the crack tip as shown in Figure 4a. The local crack tip system defines 
the tangential ( x ' or mode II) and normal (y' or mode I) coordinate directions at the crack 
tip in the deformed configuration. The terms F’ xi , F’ yi are the forces at the crack tip at 
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nodal point i in the local x and y directions respectively. The terms u' e , v' and u' p , v[* 
are the displacements at the corresponding nodal points i and l* behind the crack tip. Note 
that Gin is identical to zero in the 2D case. The extension to 3D is straightforward, as 
shown in Figure 4b and the total energy release rate Gt is calculated from the individual 
mode components as Gt =Gi+Gii +Gni ■ 

For many years the virtual crack closure technique was used mainly by scientists in 
universities, research institutions, and government laboratories, by implementing VCCT 
in their own specialized codes or in post-processing routines in conjunction with general- 
purpose finite element codes. Recently, however, YCCT was implemented, as a standard 
analysis tool, into several commercial finite element codes such Abaqus/Standard® 1 , 
Nastran™ 2 , Marc™ and ANSYS® 3 [28-32] and, therefore, has become a more frequently 
used analysis tool. 

With increasing three-dimensional simulations, the analysis along a curved or arbitrarily 
shaped front has been the focus of various new approaches. In reference [33], Li et al. 
developed a method to use VCCT in conjunction with a stepped mesh along a 
delamination front. A modified zigzag approach to approximate a crack front with an 
arbitrary shape was proposed by Liu et al. [34], A simple and efficient algorithm to trace 
a moving delamination front with an arbitrary and changing shape was developed by Xie 
and Biggers [35]. Based on the applied algorithm, a delamination front can be defined by 
two vectors that pass through any point on the front. The normal and tangent vector for 
the local coordinate system can then be obtained based on the two delamination front 
vectors. This approach allows delamination growth to be analyzed by using stationary 
meshes and it does not require the use of meshes that are orthogonal to the delamination 
front. 

A generalized method that allows arbitrary placement of the side nodes for quadratic 
elements was developed by Naim [36]. This method included both standard elements, 
with mid-side nodes, and singularity elements, with quarter-point nodes, as special cases 
of one general equation [36]. The approach also accounted for traction- loaded cracks. 
The new derivation suggested that the proper nodal forces needed for crack closure 
calculations should be the so-called nodal edge forces, rather than the global or element 
forces from standard finite element analysis results. 

Traditionally, VCCT has been used to study delamination propagation and growth where 
the delamination advance is constrained to the original interface plane. Recently, a new 
numerical method based on VCCT was proposed by Xie et al. [37, 38], which allows the 
calculation of strain energy release rates for cracks that kink. This method allows the 
evaluation of the mode I and mode II energy release rate at the tip of a kinking crack 
using only the nodal forces and displacements near the crack tip. This method may 


Abaqus/Standard® is a product of Dassault Systemes Simulia Corp. (DSS), Providence, RI, USA 

2 Nastran™ and Marc™ are manufactured by MSC. Software Corp., Santa Ana, CA, USA. NASTRAN® is 
a registered trademark of NASA 

3 ANSYS® is a product of ANSYS, Inc., Canonsburg, PA, USA 
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become of interest for analyzing the interaction of delamination propagation in one 
interface accompanied by matrix cracking, branching and subsequent delamination 
propagation at an adjacent interface. 

An interface element that incorporates VCCT was developed by Qian and Xie [39] and 
used to study an example of dynamic crack propagation under mixed mode loading. 
Using the interface element approach, VCCT can be implemented into any open 
commercial finite element code that allows user-written subroutines. The interface 
element idea is similar to the interface element proposed earlier by Mabson et al. [40-42], 

Under certain conditions, VCCT calculations yield small negative values of the mode 
components. This problem was investigated by Valvo [43], who found that the origins of 
these physically inconsistent predictions of the standard VCCT lie in the lack of energetic 
orthogonality between the crack-tip force components used to compute the mode 
components. A revised VCCT was presented which furnished a physically consistent 
partitioning of fracture modes that eliminate the calculations of negative energy release 
rate components. This was accomplished by associating the mode I and II contributions 
with the amounts of work done in a suitably defined two-step process of closure of the 
virtually extended crack. 


Recent developments for interface crack problems 

In this section, recent developments to evaluate mode mixity using VCCT and other 
approaches such as cohesive zone modeling, beam and plate modeling, and a method 
based on M-integral are presented. 


Calculating mode-mixity using the Virtual Crack Closure Technique (VCCT) 

Raju et al., in reference [26], implemented the virtual crack closure technique (VCCT) in 
a quasi three-dimensional finite element analysis. The individual mode components and 
the total energy release rate values for various bimaterial plates with isotropic, 
orthotropic, and anisotropic materials were evaluated. Reference [26] studied the problem 
of a laminated composite with an edge delamination. A typical FE model for an interface 
crack was presented in Figure 3. The total energy release rate and individual mode 
components of the strain energy release rates for a [0/35/-35/90]s laminate with 
symmetric delaminations at the 35/90 and -35/90 interfaces were calculated and analyzed. 
For all the bimaterial problems studied, the individual modes showed Aa dependence 
while the total value of G remained well defined and independent of Aa. Raju et al. [26] 
also showed that if the materials were chosen such that the Dundurs parameter p = 0 then 
the individual as well as the total G values are independent of Aa and all have well 
defined limits. 

The independence of the total value of G was also reported by several investigators [22, 
25, 27]. Specifically, Sun and Jih [27] and Dattaguru [25] observed the nonconvergence 
of the individual mode strain energy release rates and showed analytically that if the 
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oscillatory terms are neglected, the mode I and mode II values tend to be equal to one- 
half of the value of Gt i.e. Gi =Gn = Gj/ 2 at the crack tip. 

In a recent overview paper Agrawal and Karlsson [44] carefully reviewed and extended 
methods that characterize fracture at a bi-material interface, with a particular focus on 
evaluating the mode mixity using the virtual crack closure technique (VCCT). From their 
review, Agrawal and Karlsson concluded that: 

• mode mixity based on the stress intensity factor can be calculated from the 
oscillating components of the strain energy release rate obtained by VCCT, using 
the derived expressions in reference [44], 

• Sun and Jih [45], Chow and Atluri [46], and Toya [47] derived expressions for 
complex stress intensity factors using different methods. However, all these 
expressions are essentially identical and all the methods can be derived from one 
another, 

• an approach developed by Bjerken and Persson [48] yields acceptable values of 
mode mixity with significantly less computational efforts when compared to the 
other methods presented in references [46, 47] and [45], and 

• expressions suggested by Beuth [49] can also be related to the other approaches 
discussed. 

Beuth presented a mode mixity based on strain energy release rate that is da- 
independent by introducing a normalizing length parameter. However, the da- 
independent mode mixity appears to be a mathematical quantity with no 
physical meaning. Nevertheless, it is suggested that it can be a useful 
parameter if care is taken in its interpretation [44], Unfortunately, no possible 
interpretation of the quantity calculated by Beuth’ s expressions is suggested. 
Although a comprehensive summary of the mathematical relationship between the 
methods is given, no relationship of the calculated values with the commonly used 
critical values and measured mixed-mode fracture toughnesses of a material or a 
bimaterial system is suggested or provided. 


Cohesive zone modeling 

The cohesive zone modeling approach has become a widely used tool for simulating 
delaminations due to the computational convenience and ease of implementation [50-57]. 
In this approach it is assumed that a narrow zone of vanishing thickness called the 
cohesive zone exists ahead of a crack tip or delamination front. The cohesive zone 
represents the fracture process zone. The upper and lower surfaces of the narrow zone are 
held together by forces called cohesive tractions. These tractions follow a cohesive 
constitutive law (traction-separation law) that relates the cohesive tractions to the 
separation displacements of the cohesive surfaces. Delamination onset or crack growth 
occurs when the separation at the end of the cohesive zone, which represents the physical 
crack tip, reaches a critical value at which the tractions vanish. The failure process is 
controlled by displacements and stresses, which are consistent with the usual strength of 
materials theory. Thus the problem of crack tip stress singularity found in the classical 
linear elastic fracture mechanics is avoided and the singularity is effectively buried in the 
element since the crack tip is not explicitly modeled. Special finite elements, called 
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decohesion elements, with initially zero thickness, containing the traction-separation law 
can be formulated. These decohesion elements can be conveniently placed at any 
interface location in the intact model where delamination or cracking may occur. The 
modeling of an existing crack or flaw as required for VCCT to model the crack tip is not 
necessary in this approach. An extension to include modeling of fatigue crack growth 
under cyclic loading was recently presented by Naghipour et al. [58]. 

A new two-dimensional cohesive zone model, which is suitable to simulate the 
constitutive behavior of interface cracks in bimaterials, was developed by Freed and 
Banks-Sills [59]. Based upon experimental observations, a model was developed where 
the interface fracture toughness is not constant, but depends upon the mode mixity. For 
the presented approach the cohesive model may be described as an infinite set of 
cohesive laws, one for each discrete value of the phase angle yj. For each cohesive model, 
the cohesive energy and cohesive strength are adjusted with respect to the phase angle y/ 
to obtain a specific cohesive law for the current phase angle, so that only these two 
parameters are essential for this model. With this choice of parameters, the energy for 
separation Gi c is consistent with the commonly used interface fracture toughness 
criterion. 

Jin and Sun [60] reported issues, however, when cohesive zone modeling is used to 
simulate fracture at a bimaterial interface. It was found that for mixed-mode fracture in 
general, a selected single cohesive zone length may not meet the requirements that both 
the tensile and the shear stress singularities at the cohesive zone tip are cancelled 
simultaneously. The energy dissipation at the cohesive zone tip may therefore only be 
neglected when the initial stiffnesses of the cohesive model for both opening and shear 
modes are appropriately high. Also the length of the selected cohesive zone has to be set 
much greater than the characteristic length of the cohesive zone model in both cases of 
tension and shear. 


Methods based on beam and plate theory 

Recently, there appears to have been a renewed interest in analytical solutions for the 
calculation of mixed-mode energy release rates. The motivation arises from the desire to 
reduced modeling and analysis time of complicated three-dimensional finite element 
analysis. Analytical solutions based on beam and plate theory therefore become 
appealing. Bruno and Greco [61] proposed an analysis using an improved laminated plate 
model. The analysis was carried out by modeling the delaminated plate using either 
Kirchhoff or Mindlin plate assumptions for the individual layers. Comparisons with 
fracture mechanics results showed the validity of the proposed mechanical model to 
predict mode partition. In the approach suggested by Luo and Tong [62] the energy 
release rates in a cracked laminate are obtained analytically using a global - local 
method. Szekrenyes [63] proposed another analytical plate theory approach and used a 
flexible joint model to determine the displacement and stress fields in symmetrically 
delaminated, layered composite plates subject to bending. The energy release rates were 
calculated using the 3D ./-integral. 
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A methodology based on a conservative M-integral 

Banks-Sills, with several co-workers, focused on computational fracture mechanics using 
finite element analysis to study bimaterial interface fracture and compute complex stress 
intensity factors [64-70]. Over the past decade their research has also included the study 
of interlaminar fracture in composite laminates. An approach for calculating the energy 
release rate and the phase angles along the delamination front as well as measuring the 
interface fracture properties of a composite materials [68-70] was developed. 

In their analysis approach, stress intensity factors along the delamination front are 
calculated first, using a conservative M-integral formulation. The M-integral formulation 
was first presented by Yau et al. [71]. Later this formulation was extended by Yau and 
Wang [72] for interface cracks between two homogeneous isotropic two-dimensional 
bodies. More recently, a three-dimensional M-integral using asymptotic stress was 
presented by Freed and Bank-Sills [73]. 


The first term of the asymptotic stress and displacement fields, which may be expressed 
by an oscillatory singularity and a square root singularity, were obtained by means of 
Stroh and Lekhnitskii formalism s for anisotropic material. The development of the 
approach was presented in great mathematical detail in reference [73]. The stress 
intensity factors obtained from the M- integral are then used to calculate the total energy 
release rate Gt and the corresponding phase angles if/ and 0 along the front, where 


and 


Gt = 7T (K? + K?) + — K, 
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x p = arctan 
0 = arctan 
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specifically for a 0°/90° interface [74], The parameters Hi and H 2 are functions of the 
material properties. The total energy release rate and the corresponding phase angles 
from equations (10-12) are compared to the mixed-mode fracture criterion to determine 
delamination propagation. 

In general, a mixed-mode fracture criterion for a three-dimensional problem may be 
described as a relation between the three stress intensity factors or the critical energy 
release rate G c which depends on a function of the phase angles if/ and 0. It should be 
noted that the phase angle xp depends on an arbitrary length parameter L while 0 is length 
independent. The criterion used by Banks-Sills et al. [74] is based on the energy release 
rate given by 
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G c = G lc (1 + tan 2 x/j ) (1 + tan 2 0) (13) 

where G c is the critical interface energy release rate and G lc is the critical mode 1 energy 
release rate which is obtained as the average value from all tests as 


Gic 


[Re(KL iY )] 2 

K 


(14) 


Note that in equations (13) and (14) G lc is not the traditional G Ic . 

A Brazilian disk specimens, shown in Figure 5, is used to obtain the fracture parameters 
in equations (13) and (14). Note that the Brazilian disk specimens is different from the 
ASTM/ISO standard Double Cantilever Beam (DCB), End-Notched Flexure (ENF) and 
Mixed-Mode Bending (MMB) tests commonly used by the composites community. The 
fracture parameters are obtained using finite element analysis, which again is different 
from the current standards for which compliance calibration, the area method or simple 
equations based on beam theory are used to reduce the data. To date, validation of the 
method has only be established for 0°/90° and +45°/-45° interfaces. Details are given in 
reference [69]. Application examples for an embedded delamination in a cross-ply 
composite and through the width delaminations are presented in references [74, 75]. The 
validity of the method for general laminates is yet to be established. 


Additional methods to obtain satisfactory mode separation 

In this section, additional methods to obtain satisfactory individual modes are presented. 
These are the resin interlayer method, crack-tip element approach, and crack surface 
displacement extrapolation method are evaluated for their ability to obtain satisfactory 
mode separation. Each of these methods is presented along with application examples 
that demonstrate the effectiveness of the method. 


Resin interlayer method 

As an approach to eliminate the oscillatory crack-tip singularity at the bimaterial interface 
for a composite laminate, Raju et al. [26] performed analyses for a model where the 
delamination was located in a thin homogeneous isotropic layer similar to a resin layer 
that existed between the -35/90 plies of their edge-delaminated specimen. For this case 
the oscillatory component of the singularity vanished since the delamination was now 
located in a homogeneous isotropic material. The VCCT was used to calculate the total 
energy release rate GV as well as the mode I and mode II components. Results obtained 
from the model with resin layer showed that the individual mode components as well as 
the total strain energy release rates remained practically unchanged as the ratio, A a/h, of 
crack tip element length, A a, to ply thickness, h, decreased. However, in reference [76], 
no convergence for the individual mode components was observed when the interface 
crack configurations have extremely small resin layer thicknesses along with large 
stiffness differences between the layers on either side of the interface. 
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Method based on choosing A a larger than oscillatory zone 

Raju et al. [26] also reported that the individual mode components obtained from the 
interface model with the resin layer agreed well with those obtained from the models 
without resin layer when values of A a/h were either 0.25 or 0.5. The maximum difference 
between the two approaches was reported to be less than 4% for the mode I and 2% for 
the mode II components. The results suggested that the size of the delamination tip 
elements, A a, to be one-quarter to one-half of the ply thickness, h. 

Sun and Jih [27] also found that the mode components depend on the choice of the crack 
tip element length A a. However, in their analysis the mode mixity agreed well with the 
analytical solutions if the assumed crack extension was larger than the region of violent 
stress oscillation. For a calculated bimaterial constant y (also often referred to as s) for a 
AS4 composite (0.004, 0.018, and 0.038 for 0/15, 0/45, and 0/90 interfaces respectively) 
it was found that the stress-oscillating zone is negligibly small. In particular for a 90/30 
interface it was shown that the stress will not start oscillating until the distance is at the 
order of 10’ of specimen thickness [77]. Based on these results it was suggested that for 
a range of A a/a ratios (where the mode I and mode II components were found to be 
almost constant) these constant values of Gy and Gn may be used for practical purposes. 

Krueger summarized these and similar results in a comprehensive review of the state- 
of-the-art of YCCT [18]. After a careful study of the published results Krueger concluded 
that for certain values of Aa/h the variation in the individual-mode strain energy release 
rates is small as shown, in Figure 6a, and a recommended range of Aa/h was proposed. A 
value of Aa/h >1 requires smearing ply properties for plies of different orientations as 
shown in Figure 6b, which is not acceptable if individual plies need to be modeled. At 
the other extreme, a value of Aa/h < 1/20, violates the assumption of a homogeneous 
continuum since this value borders on the region where micromechanics should be 
used. As such, the recommended element size range is 1/20 < Aa/h < 1. Variations in 
mode mixity between these bounds are typically very small and should prove acceptable 
for practical applications [18]. Furthermore, the literature shows that most analysts in 
tend to use element sizes in this range. 


Application Examples 

In this subsection, two examples that illustrate the effectiveness of the methods based on 
choice of Aa larger than the oscillatory zone are presented. 

Example of a stiffened shear panel with local three-dimensional modeling^] . 

A square (1016 mm x 1016 mm) panel made of IM7/8552 carbon/epoxy tape was 
reinforced with three stiffeners made of IM7/8552 carbon/epoxy plain weave 
fabric. During manufacturing, an artificial defect of about 82 mm in length was 
introduced at the termination of the center stiffener, thus creating a delamination 
at a bimaterial interface between a -45° unidirectional tape layer and a 0° fabric 
layer. The stiffened panel was bolted to a steel picture frame (see Figure 7a) and 
subjected to shear loading. The shear loading caused the panel to buckle as shown 
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in Figure 7a. The resulting out-of-plane deformation caused skin/stiffener 
separation to initiate at the location of the artificial defect. A typical finite element 
model of this panel is shown in Figure 7b. A small section of the stiffener foot 
and the panel skin in the vicinity of the embedded defect were modeled with a 
local 3D solid model as shown in the enlargement in Figure 7b. The mixed-mode 
strain energy release rates across the width of the stiffener foot were calculated 
using VCCT. A failure index Gt/G c was calculated by correlating the computed 
total energy release rate Gt and mode mixity Gn/Gj with the mixed-mode failure 
criterion of the graphite/epoxy material similar to the one shown in Figure 1. 

The objective of the research was to study the effect of the fidelity of the 
local 3D finite element model on the computed mixed-mode strain energy release 
rates and the failure index [78]. This study also included the introduction and 
modeling of a resin layer at the delaminated interface between the top skin ply (-45° 
unidirectional tape) and the bottom of the flange (0° fabric layer). As the 
delamination exists in a homogeneous isotropic material (the resin interface layer) 
the oscillatory singularity at the crack tip of the bimaterial - as mentioned above - 
does not exist. For convenience, the delamination was placed centrally within the 
resin layer. Therefore, the resin was modeled with two layers of elements with the 
delamination located between the two layers. Further details are discussed in 
reference [78]. The failure index Gt/G c distribution calculated along the 
delamination front, s, was compared to the distribution calculated from models 
without a resin layer as shown in Figure 8. For the models without a resin layer, a 
crack tip element length, A a, had been chosen (approximately one ply thickness, 
Aa/h=l ) in the range over which the mode mixity exhibits a reduced sensitivity to 
the value of A a as discussed above. The results from both approaches were in good 
agreement with each other and showed that the approaches yielded nearly equivalent 
results, confirming the results from Raju et al. [26] for the edge delamination 
problem. Based on the comparison, local models without the resin layer were used 
for all subsequent studies to reduce model complexity and analysis time [78]. 

Example of a multi-directional Single Leg Bending (SLB) specimen 
To confirm the assumptions made about the range over which the mode mixity 
exhibits a reduced sensitivity to the value of A a, as discussed above, a Single Leg 
Bending (SLB) specimen as shown in Figure 9a was studied in reference [18]. 
The SLB specimen was designed for mixed-mode fracture toughness testing and 
exhibits a nearly constant mode ratio of Gii/Gt=0A (40% mode II)[79]. A 
specimen with multidirectional layup was analyzed, where the delaminated 
interface was located between a +30° and -30° ply [80, 81]. The three- 
dimensional model of the SLB specimen is shown in Figure 9b. Along the length 
of the model, a refined mesh of length c was used in the vicinity of the 
delamination front. The influence of mesh size on computed mixed mode strain 
energy release rates was studied by keeping the length of the refined zone, c, 
constant and increasing the number of elements, n, in this zone, as shown in the 
detail of Figure 9b. The corresponding values of relative element size, A a/h, and 
relative crack closure length, A a/a, are shown in the insert. 
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The influence of mesh refinement on the total strain energy release rate 
distribution Gt across the normalized width, y/B, is shown in Figure 10. Only very 
long elements {Aa=c/n=\ mm; Aa/h=S ) need to be avoided since the model is too 
coarse to yield any sensible results. The results for the other models basically 
agree very well with one another. The distribution of the mixed mode ratio Gh/Gt 
is fairly constant at about 35% mode II across almost the entire width of the 
specimen as shown in Figure 11. For the range studied, only a small dependence 
of computed mixed mode ratio on element size Aa/h was observed. 

Similar results were also obtained in a related investigation, in which mixed-mode 
energy release rates were computed for an ENF specimen with the same layup 
and interface location as the SLB specimen [80, 82], This study indicated that the 
computed energy release rates, Gu and Gm, did not exhibit a significant variation 
when Aa/h was kept within the suggested range. The results confirm the 
suggestion made above that element lengths be chosen within lower and upper 
bounds. 

These methods based on choice of Aa larger than the oscillatory zone are a convenient 
practical engineering approach based on experience. However, this simple approach may 
be limited to cases where the material mismatch is insignificant and the bimaterial 
constant, y, is small. For general interface cracks along a bimaterial interface, where there 
is a significant mismatch in material properties and stiffnesses, this issue is much more 
difficult to resolve. These significant mismatches are observed for coatings on substrates 
in the electrical chip industry [8] or for sandwich structures with stiff face sheets and 
compliant cores [76]. 


Crack tip element (CTE) approach 

The most comprehensive of all methodologies based on plate theory is the so-called 
crack tip element (CTE) approach developed by Davidson with several co-workers [83- 
87] over the past 20 years. His overall prediction methodology includes the analysis of 
the energy release rate and the mode mixity along the delamination front in a structure as 
well as the development of a suitable mixed-mode fracture criterion [88-90]. The 
proposed mixed-mode fracture criterion is based on fracture toughness testing using 
simple standard type specimens such as the DCB, ENF and MBB specimen. The 
proposed mixed-mode fracture criterion is similar to the criterion shown in Figure 1. The 
Single Leg Bending (SLB) specimen, for mixed-mode testing at approximately constant 
40% mode II was also introduced [79]. The methodology provides a straightforward 
practical solution for engineering applications. 

A linear elastic crack-tip element for two-dimensional analysis was developed first 
and then extended and refined for practical applications [83] [84], Analytical expressions 
were derived for the total energy release rate and mode mixity based on plate theory force 
and moment resultants near the crack tip as shown in Figure 12a. The element may be 
used for cracks within or between homogeneous, isotropic, or orthotropic layers, or for 
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delaminations in laminated composites. After successful application to two-dimensional 
problems, the method was extended and a 3D crack tip element was developed as shown 
in Figure 12b [85, 86]. Classical plate theory was used to derive the equations for the 
total energy release rate and mode mix. An additional mode mixity parameter, Q, is 
required to complete the mode mixity decomposition. This parameter depends upon the 
elastic and geometrical properties of the materials above and below the crack plane, but 
not on the loading. A relatively simple finite element model may be used to determine the 
mode mixity parameter. Pre-calculated values of Q were also presented for a large 
number of cases [83]. For bimaterial interfaces where an oscillatory singularity occurs at 
the crack tip, an approach was also described which allows a unique, physically 
meaningful value of mixed mode ratio to be defined [83]. Davidson refers to this 
approach as the / 3=0 approach, where the generalized Dundurs parameter [1 is taken equal 
to zero only in the finite element model that is used to to calculate Q. Setting [1=0 causes 
the bimaterial constant y in equation (1) (also often referred to as s) to be zero and 
eliminates the oscillating singularity at the crack tip [91, 92]. Setting [1=0 is accomplished 
by altering the Poission’s ratios v 13 or v 23 of one of the layers at the interface. The elastic 
moduli, however, need to be altered in a manner such that [1=0 while the total energy 
release and other significant parameters remain unaffected. With [1=0 and y =0 the mode 
mixity parameter Q can be calculated as 


sin ft 


-V5 

WcVCii 


(15) 


Here, N c is the concentrated crack tip force and cn, is a function of the material 
properties and layup of the plate. The procedure is described in full detail in reference 
[83], 


For practical purposes, the energy release rate distribution along the delamination front is 
obtained from each crack tip element along the delamination front. Each crack tip 
element is made of a generic set of four elements as shown in Figure 12c. The total 
energy release rate Gt is obtained by a plate-theory-based crack closure procedure as 

2 

G t =- ^(AtyA^+AMiAic?);, i = 1,2,6 (16) 

7 = 1 

In order to obtain the quantities A A, , A M t , A e° and A K t for each crack tip element along 
the front, the plate theory forces, {N\,N%N( } \ and moments , obtained at the 

centroids of the four elements are first used with laminated plate theory to obtain the 
strains, e, and curvatures, k, in each of the cracked regions. The mode I and mode II 
components can be calculated as 

\ 2 

Gi= - [-Jc[[N c sin D. + Jc^M c cos(D. + Y)] (17) 
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where 


( 18 ) 


\ 2 

G n = - l-Jc^Nc cos H + Jc^M c sin(H + T)] 

C12 

sin T = : (19) 

V C 11 C 22 

Here, V c and M c are the concentrated crack tip force and moment, respectively. The c\\, 
C 22 and C 12 terms are functions of the material properties and layup of the plate which are 
obtained from the inverse of the reduced stiffness matrix used in classical laminated plate 
theory [85, 86, 88, 93]. The mode III component is calculated from these results as 

Gm= G — Gj — G u (20) 

Additional details on the 3D CTE analyses are presented in reference [85]. In a detailed 
application example, delamination onset from an embedded defect is predicted using the 
3D CTE analysis [86]. The defect is located between the stiffener foot and the skin of a 
hat stiffened specimen subjected to tension and bending. Good correlation with 
delamination growth measured in related experiments was observed. 

Recently, Mikulik et al. [94] investigated the accuracy of VCCT and CTE to predict the 
initiation and growth of delaminations for the development of a specimen with multiple 
delaminations at several interfaces. Their investigation may be regarded as an 
independent verification of the CTE approach. The predictions correlated very well with 
the test data and both numerical approaches demonstrated the capability to identify the 
critical delamination interface and accurately capture the initiation load levels. Mikulik et 
al. concluded that fracture mechanics based approaches can be used for the prediction of 
multi-level delamination propagation. 


Crack Surface Displacement Extrapolation (CSDE) method 

Sandwich constructions used in aerospace vehicle structures generally consist of a low 
density core between two stiff surface layers (face sheets) to achieve lightweight 
structures with high stiffness and strengths. Flaws in the form of disbonds between the 
face sheet and the core may reduce the load-bearing capacity of the structure because of 
the loss of load transfer between the face sheet and the core. Such a disbond generally 
constitutes an interfacial crack between two widely dissimilar materials [95]. Smith and 
Shivakumar [96] used the total energy release rate Gt to analyze their Cracked Sandwich 
Beam specimen to avoid issues with respect to mode separation that could occur when 
using VCCT. Wang and Zhang [97] presented a new analytical solution for the energy 
release rate and the phase angle at an interface crack in a sandwich based on interface 
fracture mechanics solutions and composite beam theory. 

A promising method for the prediction of disbond propagation at the sandwich face 
sheet/core interface was suggested by Berggreen et al. [98, 99]. The CSDE method is 
based on the complex stress intensity factor K defined in equation (5) and the so-called 
Crack Surface Displacement technique proposed by Matos et al. [100]. 

The total energy release rate Gt and the mode mixity phase angle i//y proposed in [92] 
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x/j k = arctan 


Im (Kh iy ) 
R e(Kh}Y') 


( 21 ) 


can be expressed in terms of the relative crack flank displacements S x and d y shown in 
Figure 13a as 


where h is the characteristic length of the problem and is chosen here as the face 
thickness. The constants Hn, H22, and the oscillatory index y are dependent on the 
materials at the interface. 

The mode-mixity and energy release rate can be calculated using equations (22, 23). The 
relative nodal displacements of the crack flanks are obtained from a finite element 
analysis. Therefore, both the mode-mixity, y/K, and Gt will be affected by the oscillations 
near the crack tip. However, these effects are believed to be small away from crack tip 
zone and therefore have no practical importance. This argument is similar to the 
assumption made earlier for a practical solution of the mode mixity which can be 
obtained if the crack tip element length A a was larger than the region of violent stress 
oscillation [27]. As shown in Figure 13b, where the total energy release rate Gt (solid 
blue line) was plotted for crack flank node pairs (solid blue circles), the oscillations 
vanish with increasing distance from the crack tip. A similar graph can be plotted for the 
mode-mixity. To avoid the oscillations in the numerical error zone, only the relative 
crack flank displacements obtained in the evaluation zone are used in the Crack Surface 
Displacement Extrapolation method. The results obtained are simply linearly extrapolated 
into the crack tip (red star in Figure 13b). The extrapolation can be performed separately 
for the total energy release rate, Gt, and the mode-mix, iff K , but normally it is sufficient to 
do this numerical extrapolation once for Gt and then reuse the linear extrapolation 
parameters to obtain the mode mixity curve. Results for a rectangular sandwich plate with 
a circular disbond were presented where the computed energy release rate and mode 
mixity for the structure were compared to measured fracture toughness values [98, 99]. 

In another example where the fracture of a sandwich specimen loaded with axial forces 
and bending moments was analyzed, the CSDE method was used to calculate the phase 
angle [101, 102], Face sheet/core fatigue crack growth in foam-cored sandwich 
composites was examined by Manca et al.al [103] using the mixed mode bending (MMB) 
test method. The phase angle, y/, and energy release rate, G, were both determined from 



( 22 ) 


and 



(23) 
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two-dimensional plane strain finite element results using the CSDE method. The method 
was also applied in a recently proposed accelerated simulation scheme for fatigue crack 
growth in a bimaterial interface, which was demonstrated to simulate fatigue crack 
growth in the face sheet/core interface of a sandwich beam [104], 


Discussion 

Although many new methods and approaches are presented in the literature for the 
interface crack problem, only a few approaches appear to be of practical use and capable 
of application to general cases and real world structural problems. Until now, the 
community at large has not yet found a consensus on an approach or a set of approaches. 
New methodologies need to provide analysis tools to determine the mixed-mode fracture 
parameters along an arbitrarily shaped crack or delamination front located at a bimaterial 
interface in a three-dimensional solid. In order to be successful, these new methodologies 
simultaneously need to recommend methods that are suitable for fracture toughness 
testing. Relatively simple specimens and uncomplicated test procedures should be 
developed to generate the data required for use in a mixed-mode fracture criterion. Any 
mixed-mode fracture criterion developed can then be used in conjunction with the 
analysis results to determine crack onset and growth. 


Concluding Remarks 

Cracks at bimaterial interfaces are common in composites, microelectronic materials, 
film interfaces, etc. To be able to predict the fracture behavior and perform damage 
tolerance evaluations of these materials and structural components, a thorough 
understanding of behavior of cracks located at bimaterial interfaces is needed. Williams 
showed more than a half a century ago, that in addition to the classical square root 
singularity at the crack tip, there exists an oscillatory singularity for cracks located at a 
bimaterial interface. The oscillatory singularity leads to some unsatisfactory 
consequences. Several investigators over the past three decades showed that when 
numerical methods, such as the finite element method, are used to evaluate the total and 
individual mode strain energy release rates, the individual modes do not show 
convergence as the mesh size is refined near the crack tip. The intent of the review was 
to renew the awareness of the oscillatory singularity at the crack tip of a bimaterial 
interface and the problems that occur when calculating mode mixity using numerical 
methods such as the finite element method in conjunction with the virtual crack closure 
technique. 

Established approaches to overcome the non-convergence issue of the individual mode 
strain energy release rates were reviewed. In the recent literature many attempts to 
overcome the nonconvergence issue have been developed. Among the many approaches 
found only a few methods hold the promise of providing practical solutions. These are 
the resin interlayer method, the method that chooses the crack tip element size greater 
than the oscillation zone, the crack tip element method that is based on plate theory and 
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the crack surface displacement extrapolation method. Each of the methods was validated 
on a very limited set of simple interface crack problems. However, their utility for a wide 
range of interfacial crack problems is yet to be established. 
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. Mixed-mode failure criterion for a typical carbon/ epoxy material. 




Figure 2: Crack located at a bimaterial interface. 
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b. VCCT for eight node 3D solid elements. 

Figure 4. Virtual Crack Closure Technique (VCCT) [17,18], 



Figure 5: Sketch of a Brazilian disk specimen containing a composite 
strip with a crack along a 0/90 interface [69,74], 



Aa: Element length 
h: Ply thickness 


(a). Dependence of computed energy release rate components on element size near the crack tip. 
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(b). Upper and lower bounds for element size at crack tip 
Figure 6: Influence of element size on computed mode components . 
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(a). Buckled composite panel under shear loading with picture frame. 
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(b). Finite Element model of stiffened panel and load frame. 
Figure 7. Stiffened composite panel [78] . 
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location along delamination front, s 

8 . Computed failure index across the width of the stiffener for 
models with and without a resin layer [78]. 



{ 2 

2L 

a 


Dimensions 

25.4 mm 
2.032 mm 

2.032 mm 

177.8 mm 
34.3 mm 


[±30/0/-30/0/30/0 4 /30/0/-30/0/-30/30//-30/30/30/0/30/0/-30/0 4 /-30/0/30/0/±30] 

Tl 


Delamination 


(a). SLB specimen 
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(b). Three-dimensional finite element model of a SLB specimen with detail. 

Figure 9. Finite element analysis of a Single Leg Bending Specimen (SLB) [18, 80, 81 ]. 
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Figure 10. Influence of number of elements in refined section on computed total strain energy 
release rate distribution across the width of a SLB specimen [18]. 



Figure 1 1 . Influence of number of elements in refined section on mode ratio distribution 

across the width of a SLB specimen [18]. 



(a). Two-dimensional crack-tip element [83 ]. 
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(c). Four finite plate elements adjacent to the crack front in a double plate model [88] . 
Figure 12: Crack-tip element approach suggested by Davidson. 



(a). Sandwich bimaterial interface and crack tip geometry. 



(b). Detail of evaluation zone. 

Figure 13: Schematic of the crack surface displacement extrapolation method [98, 99]. 



